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By breaking the reflection and the axial symmetries simultaneously, we developed multi- 
dimensional constraint relativistic mean field (MDC-RMF) models. In these models, the nuclear 
shape is assumed to be invariant under the reversion of x and y axes, i.e., the intrinsic symmetry 
group is V4, and all shape degrees of freedom Px^ with even fj, (/320, /322, Pso, /J32, /340, ■■•) are 
included self-consistently. The RMF functional can be one of the following four forms: the meson 
exchange or point-coupling nucleon interactions combined with the non-linear or density-dependent 
couplings. The pairing effects are taken into account with the BCS approach. In this paper the 
formalism for MDC-RMF models is presented in details. Potential energy surface of ^''"Pu is il- 
lustrated for numerical checks and for the study of the effects of triaxiality on the fission barriers. 
Potential energy curves of actinide nuclei around the first and second fission barriers are studied 
systematically with emphasis on the second ones which are lowered considerably by the triaxial 
deformation. We conclude that it is important to include the reflection asymmetric and non-axial 
shapes simultaneously for the study of potential energy surfaces and fission barriers of actinide nuclei 
and of those in unknown mass regions such as superheavy nuclei. 

PACS numbers: 21.60.Jz, 24.75.+i, 25.85.-w, 27.90.-fb 
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I. INTRODUCTION 

The occurrence of spontaneous symmetry breaking 
leads to nuclear shapes with a variety of symmetries [1, 2]. 
A lot of nuclear phenomena are connected with the nu- 
clear deformation, including the small and large ampli- 
tude collective motions, e.g., the rotation and the fis- 
sion [3, 4]. There are also many unresolved problems 
concerning the motion or change along the nuclear shape 
degrees of freedom. For example, in recent years lots of 
efforts were devoted to explore nuclear shape phase tran- 
sitions [5-8]; one way to investigate shape phase transi- 
tions is to examine nuclear potential energy surfaces and 
their changes with the nucleon number [9-13]. 

The shape of a nucleus can be described by the 
parametrization of the nuclear surface or the nucleon den- 
sity distribution. There are mainly two kinds of ways to 
parametrize the nuclear shape [14, 15]. One of them is 
the two-center or two-center-like parametrization [16-20]. 
The other way is of one-center and to make a multipole 
expansion. 
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where /3a/^'s are deformation parameters. The way of 
multipole expansion is usually used in mean field cal- 
culations. The majority of observed nuclear shapes is 
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of spheroidal form which can be described by the axial- 
quadrupole deformation parameter /32o though in early 
years the Nilsson pcrturbed-spheroid parameter €2 was 
usually adopted for numerical convenience [21]. The non- 
axial-quadrupole (triaxial) deformation /322 (or, equiva- 
lently, 7) has been predicted in atomic nuclei for quite 
a long time and it may play important roles in super- 
heavy nuclei [22]. A static triaxial shape manifests it- 
self by the wobbling motion and chiral doublet bands; 
both have been extensively studied from experimental 
and theoretical sides [23-27]. Particularly, the existence 
of multiple chiral doublet (M^D) bands in one nucleus 
was predicted [28]. Recently two distinct chiral doublet 
bands were observed in ^-^^Ce [29] which for the first time 
provides a convincing verification of the prediction. The 
softness with respect to the gamma distortion also results 
in many interesting nuclear phenomena [30-35]. 

The octupole shapes with A = 3 are predicted to ex- 
ist in nuclei in several mass regions [36]. The low- lying 
negative parity bands observed in actinidcs and some 
rare-earth nuclei are related to the reflection asymmetric 
shapes [37-43]. It was also found that the octupole de- 
formation is responsible for the Ra puzzle [44], i.e., the 
anomalous enhancement of the residual proton-neutron 
interactions for Ra isotopes around N = 135 [45, 46]. 
The influences of the non-axial octupole /332 deforma- 
tion on the low-lying spectra have been also traced out 
theoretically and experimentally [47-56]. For example, 
reflection asymmetric shell model calculations revealed 
that the observed low-energy 2~ bands in N = 150 nu- 
clei [57] are caused by the P32 deformation [58]. Indeed, 
strong 132-correlations were found in some N = 150 iso- 
tones from multi-dimensional constraint covariant den- 



sity functional theories (MDC-CDFT) [59]. Experimen- 
tal efforts have also been made to explore the tetrahedral 
symmetry in nuclei [60, 61]. 

Deformations of higher multipoles with A > 3 are im- 
portant to different extents. The hexadecapole deforma- 
tion, /340 or £4, has been included in deformed mean field 
potentials since 1960s, see, e.g., Ref. [62]. The important 
effects of the high order deformation /Sgo or eg on the 
angular momentum alignments and dynamic moments of 
inertia in superheavy nuclei are also discussed [63, 64]. 

The shape degrees of freedom are important not only 
for the ground states or small amplitude collective mo- 
tions, but also for large amplitude collective motions such 
as fission. Since the discovery of the nuclear fission [65], 
the description of the fission process has been a difficult 
and challenging task. The fission dynamics are mostly 
determined by the barriers which prohibit the dissolving 
of the nucleus [66] . In order to study the fission problem, 
one should have very accurate information about the fis- 
sion barrier, i.e., the height and width, or more precisely, 
the shape of the fission barrier [67, 68]. Particularly, 
to explore the island of stability of superheavy nuclei 
(SHN) [69-72], it is more and more desirable to have ac- 
curate predictions of fission barriers of SHN [73-81]. Un- 
til now many popular nuclear structure models have been 
employed to calculate the nuclear fission barriers, includ- 
ing the macroscopic-microscopic (MM) model [18, 73, 82- 
86] , the extended Thomas- Fermi plus Strutinsky integral 
(ETFSI) method [87], the Hartree-Fock or Hartree-Fock- 
Bogoliubov method with the Skyrme force [88-93] and 
the Gogny force [94, 95], and the covariant density func- 
tion theory [88, 96-102]. 

Besides /32o which describes the elongation of a fissile 
nucleus and /34o which is relevant to the size of a neck, 
many other shape degrees of freedom are also crucial for 
determining the shape of fission barriers and the fission 
path. Let's take actinides as examples. Due to the shell 
effects, actinide nuclei are characterized by a two-humped 
fission barrier [14] and there may also exist a third min- 
imum [86, 103-107]. It has long been known from the 
MM model calculations that the inner fission barrier is 
lowered by the non-axial-quadrupole deformation [108- 
110] and the outer barrier by the reflection asymmetric 
(RA) shape [111]. Later, the important roles played by 
the non-axial-quadrupole deformation and octupole de- 
formation were confirmed in the non-relativistic [112] and 
relativistic [113] density functional calculations, respec- 
tively. Therefore what is usually done is to consider the 
triaxial but reflection symmetric (RS) shapes for the in- 
ner barrier and axially symmetric but RA shapes for the 
outer one [73, 94, 114]. The non-axial octupole deforma- 
tions are considered in both the MM models [115] and 
the non-relativistic Hartree-Fock theories [116]. 

In recent years, the nuclear CDFT has been very suc- 
cessful in the description of both the ground states and 
the exciting states of the nuclei ranging from light to 
superheavy regions [117-124]. Particularly, the CDFTs 
have been used to study the potential energy surfaces 



(PES) and the fission barriers of heavy and superheavy 
nuclei [100, 101, 120, 125-128]. We have developed MDC- 
CDFTs by breaking the reflection and the axial symme- 
tries simultaneously. In these theories, all shape degrees 
of freedom Px^ with even ^, e.g., /32o, /322, Pso, Ps.2, 
/340, • • • , are included self-consistently. The functional 
can be one of the following four forms: the meson ex- 
change or point-coupling nucleon interactions combined 
with the non-linear or density-dependent couplings. For 
the particle-particle channel, either the BCS approach 
or the Bogoliubov transformation has been implemented 
and the pairing force could be 5-force or the finite-range 
separable force [129-131]. For convenience, the MDC- 
CDFT with the BCS approach for the pairing is named as 
MDC-RMF model and that with the Bogoliubov trans- 
formation as MDC-RHB one. Due to the heavy com- 
putational burden, MDC-RHB theories have not been 
used to do multi-dimensional constraint calculations yet 
(here "multi" means three or more). In Ref. [125, 126] 
we calculated the potential energy surfaces of actinide 
nuclei in the (/320j 1^22, Pao) deformation space with the 
MDC-RMF models and we found that the triaxiality also 
plays an important role upon the second fission barriers. 
In this paper, we will present the detailed formulae for 
MDC-RMF models and some results of actinide nuclei. 

The paper is organized as follows. The formalism of 
the MDC-RMF models wiU be given in Sec. II. In Sec. Ill 
we present the numerical details and results on the PES's 
of the actinide nuclei. A summary is given in Sec. IV. 



II. FORMALISM OF MULTI-DIMENSIONAL 

CONSTRAINED RELATIVISTIC MEAN FIELD 

MODELS 



In the CDFTs a nucleus is treated as a composite of 
nucleons interacting through exchanges of mesons and 
photons. The effects of mesons are described either by 
the mean fields or by the point-like interactions between 
the nucleons. Meanwhile, the density dependence of the 
coupling constants [132, 133] or the non-linear coupling 
terms [134-136] are introduced to give correct satura- 
tion properties of the nuclear matter. Accordingly, the 
form of the covariant density functional can be one of the 
following four: the meson exchange or point-coupling nu- 
cleon interactions combined with the nonlinear or density 
dependent couplings. Most of the computational efforts 
are devoted to the solving of the Dirac equation and the 
calculation of the various densities which are common for 
all these CDFTs. In this section, we mainly present the 
formalism of the CDFT with the non-linear point cou- 
plings (NL-PC). For more details of CDFTs, the readers 
are referred to Refs. [117-124]. 

The starting point of the covariant density functional 
with NL-PC is the following Lagrangian: 
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are the hnear couphng, non-hnear coupling, derivative 
coupling, and the Coulomb part, respectively. M is the 
nucleon mass, as, ay, ars, cerv, Ps, Is, Iv, 5s, Sy, 
Sts, and Stv are coupling constants for different chan- 
nels and e is the electric charge, ps, Pts, jv, and jtv 
are the iso-scalar density, iso-vector density, iso-scalar 
current, and iso-vector current, respectively. The various 
densities and currents are defined as: 
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(4) 



Starting from the Lagrangian, using the Slater deter- 
minants as trial wave functions, and neglecting the Fock 
term as well as the contributions to the densities and 
currents from the Dirac sea, we can derive the relativis- 
tic mean field equations with the variational principle: 
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where 



[p-V{r)]+^[M + Sir)] + Vo{r), (6) 



is the single particle Hamiltonian and 

S = asps + OiTSPTS ■ f + PsPs 
+ Ss^ps + Sts^Pts ■ T, 
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are the scalar and vector potentials, respectively. In the 
present work we suppose that the states are invariant 
under the time-reversal operation, which means that all 
the time-odd or the vector components of the currents 
and the potentials vanish. In this case the single par- 
ticle Hamiltonian has the time-reversal symmetry which 
simplifies the calculations. For convenience two vector 
densities are defined as the time-like components of the 
4-currents jv and Jtv- 



Pv{r) = jv{r), pTv{r) = JTv{r) 



(9) 



which are the only non-vanishing components. 

It is customary to solve the deformed CDFT by ex- 
panding the axillary single particle wave functions on 



a group of complete basis, e.g., the harmonic oscilla- 
tor (HO) basis [137, 138] or the Woods-Saxon (WS) ba- 
sis [139-143]. A reflection asymmetric relativistic mean 
field (RAS-RMF) approach has developed by using a ba- 
sis in a two-centre harmonic-oscillator potential [144]. In 
our MDC-CDFTs, the single particle wave functions and 
various densities are expanded on an axially deformed 
harmonic osciUator (ADHO) basis [137, 138, 145]. The 
ADHO basis are defined as the eigen solutions of the 
Schrodinger equation. 



2M 



y^ + VB{z,p) 



<^^{ra)^E^'^„{ra), (10) 



where 
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is the ADHO potential and Wz and ojp are the oscillator 
frequencies along and perpendicular to the symmetry z 
axis, respectively. The solution of Eq. (10) reads 
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where (j)n,{z) and R™'^{p) are the HO wave functions, 
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Xs^ is a two component spinor and Ca is a complex 
number introduced for convenience. Harmonic oscilla- 
tor lengths hz and bp are related to the frequencies by 
hz ~ Xj \JMu)z and hp = 1/ yjMup. The corresponding 
eigen energy is Ea = ujp{2np + |?ti;| + l)+ujz{nz + \) and 
the major quantum number is N^ = 2np + \mi\ -\- Uz- 

These basis are also eigen functions of the z compo- 
nent of the angular momentum, jz with eigen values 
K ^ mi + rris. For any basis ^^(rcr) the time re- 
versed basis is defined as ^ai^cr) = T^ai'i'O'), where 
T = ifjyK is the time reversal operator and K is the com- 
plex conjugation. Apparently we have Ka = —K^ and 
TTg = tTq where tt = ±1 is the parity. The deformation 
of the basis /3basis is defined through the relations uJz = 

f-y^/3basisj and lOp = wq cxp K 

where wq = {'jJz^'iY is the frequency of the correspond- 
ing spherical oscillator [137, 138]. 

These basis form a complete set for expanding any two 
component spinors. For a Dirac spinor with four compo- 
nents. 
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where the sum runs over all the possible combination 
of the quantum numbers a = {nz,nr,mi,ms}, /" and 



gf are the expansion coefficients. In practical calcu- 
lations the sum in Eq. (14) arc truncated as riz/Qz + 
{2np + \m\)/Qp < A^ph, where A^ph is an integer and 
Qz = max(l, 6z/6o) and Qp = max(l,6p/6o) are con- 
stants related to the oscillator lengths bo, bz, and bp. 
Clearly in any case the number of the basis included are 
no less than that obtained from the major shell trunca- 
tion N < Nph. 

If a nucleus is invariant under the rotation around z 
axis and the spatial reflection, the angular momentum 
projection and the parity are conserved and the RMF 
equation (5) can be decomposed into blocks characterized 
by the quantum numbers K and tt. Usually only half of 
the space with K^ > is considered due to the time- 
reversal symmetry. 

Now let us turn to the general case that the axial sym- 
metry as well as the spatial reflection symmetry arc bro- 
ken. If we still expand the spinors on the ADHO basis, 
different components with different K and tt are mixed 
together, that is, we must diagonalize a larger single par- 
ticle Hamiltonian matrix with non-zero matrix elements 
between two basis states with different K and tt. Nev- 
ertheless, even in this case we still can have one symme- 
try operator that makes the Hamiltonian matrix block- 
diagonal. Due to the axial symmetry of the basis, it is 
convenient to choose the simplex operator S = ze~"^^= 
as a conserved quantity. Note that for a fcrmion system 
with a half-integer spin, S* is a Hermitian operator and 
S^ = I. This operator corresponds to the rotation by 
TT around the z axis thus leaving the nuclear mean field 
invariant. S is also a good quantum number for the basis 
with S^a = (— 1)^^°~^^^*&Q, which means that the basis 
$a with Ka = +1/2, -3/2, 5/2, -7/2, • • • span the sub- 
space with S ~ I, while their time- reversed states span 
the one with 5* = — 1. Note that now the blocks with 
K = -hl/2 and K = -3/2 (not that of i^ = +3/2) are 
mixed. Remembering that for a nucleus with the time 
reversal symmetry only the basis with S = 1 are used in 
the expansions, for such basis we set Ca ~ I when ex- 
panding the large components and Ca — i for the small 
ones. The basis with Sa = —1 are obtained by simply 
applying T on them. Furthermore, for systems with the 
time-reversal symmetry, it is only necessary to diagonal- 
ize the matrix with S = I and the other half is obtained 
by a time reversal operation on the obtained single par- 
ticle wave functions. 

For deformed nuclei with the V4 symmetry, we expand 
the potentials V{r), S{r), and the densities in Eq. (4) in 
terms of the Fourier series. 
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By substituting this expansion into the symmetry condi- 
tions, it is easy to see that fii = fu= fp, and /« = for 



odd n. Thus the expansion (15) can be simplified as 
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where 
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are real functions of p and z. The details for calculating 
the matrix elements of the Dirac Hamiltonian and various 
densities and their derivatives are given in the Appendix. 
For open shell nuclei the pairing interaction becomes 
crucial and must be included. For example, the calcu- 
lated fission barriers depend very much on the form and 
strength of the effective pairing interactions [146]. In 
realistic calculations several methods have been used to 
treat the pairing effects, e.g., the BCS approach, the Bo- 
goliubov transformation, and the particle number con- 
serving method [147-151]. Here we show the formulas for 
the BCS approximation with a S pairing force which are 
incorporated in the MDC-RMF models. In the particle- 
particle channel, the gap equation reads [2], 
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where e^ — e^ — A and A is the Fermi energy. The total 
pairing energy is 



^pair = -^Vo I dVK*(r)K(r). 



(19) 



To obtain the PES one can perform a constraint calcu- 
lation which is equivalent to adding an external potential 
during the iteration [2] . The quadratic constraint method 
is usually used. 
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where C\^ is the spring constant and mx^'a are desired 
moments. With this method the calculation always con- 
verges to a deformation point on the PES other than the 
desired one. To overcome this shortcoming and to get a 
PES with equally distributed points, we use a modified 
linear constraint method. The Routhian reads. 



E' — -Ermf + 2_^ -^Cxf^Qxfj,, 



(21) 
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with the variables Ca^'s changing their values during the 
iteration through the following relation. 
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where /3ap is the desired deformation, k\^ is a constant, 

(n) 

and Cj^ is the value at the nth step. We found that this 
constraint method works well in our multi-dimensional 
calculations. 

The total energy of the nucleus reads 

+ 2'^spl + -j^oLvPv + ^o-tsPts + ^oitvPtv 
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where Ec.m. is the center of mass correction which can 
be calculated either from the quasi-particle vacuum. 
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where P is the total linear momentum and A is the nu- 
clear mass number or in the oscillator approximation, 

-Bern. =-| x41 X Ai/3 MeV, (25) 

depending on the effective interactions used in the RMF 
functional. 

The intrinsic multipole moments are calculated from 
the densities by 



Qxt. = j d'rpviryY^^in), 



(26) 



where Y\^{^) is the spherical harmonics. The deforma- 
tion parameter (3\^j^ is obtained from the corresponding 
multipole moment by 
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where R is the estimated nuclear radius R — 1.2x^^/'^ fm 
and N is the number of protons, neutrons, or nucleons. 



III. RESULTS AND DISCUSSIONS 
A. Numerical details 

In this section we present the numerical details and 
some illustrative calculations of MDC-RMF models. The 
potentials and densities are calculated on a spatial lat- 
tice of which mesh points in the p and z directions are 
designed in a way that the Gaussian quadrature can be 
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FIG. 1. (Color online) The potential energy curve of ^^''Pu 
calculated with different number of the oscillator shells A'ph 
from the MDC-CDFT calculations. The axial symmetry is 
imposed. The results calculated with A^ph = 16, 18, and 20 
are depicted by solid, dashed, and dotted lines, respectively. 
The four sub-figures show the detailed structure of the PES's 
near the ground state, the inner barrier, the isomeric state, 
and the outer barrier. The width of each sub- figure is 0.1 and 
the height is 1 MeV. 



made and those for the azimuthal angle are equally dis- 
tributed. Since we keep the mirror reflection symme- 
try with respect to the a; = or j/ = planes, only 
mesh points with positive x and y are considered. For 
the azimuthal angle, more than 10 mesh points for light 
nuclei and about 20 mesh points for heavy nuclei are 
used, which are enough for most of practical applications. 
In two special cases the number of the mesh points can 
be even reduced: (1) for axial symmetric nuclei the az- 
imuthal degree of freedom vanishes and (2) for reflection 
symmetric nuclei the mesh points with z < can be 
omitted. With these choices the values of the localized 
fields and potentials in the full lattice space can be simply 
obtained by symmetry transformations like the rotations 
or space reflection. Nevertheless, the Coulomb field must 
be treated carefully due to its long range nature. 

The calculated physical observables should converge as 
the truncation A'ph — >■ 00. In Fig. 1 we show the poten- 
tial energy curve of ^^°Pu calculated with different trun- 
cations, A'^ph = 16, 18, and 20. The effective interaction 
PC-PKl [152] is used. For simplicity the calculations 
are restricted to the axially symmetric case. As is found 
in most of earlier calculations, the results show a typi- 
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against the basis deformation. Furthermore, the results 
with iVph = 16, 18, and 20 coincide with each other. This 
conclusion holds also for the second barrier, /32o ~ 1-3, 
except that two points with small basis deformations and 
A^ph = 16 are very high. However, if we choose properly 
the basis deformation, 16 shells are enough for calcula- 
tions around the barriers. For even larger deformation, 
/320 = 2.0 or 3.0; the results for A^ph = 16 deviate from 
that for A^ph = 18 and 20; the differences between the 
results for TVph = 18 and 20 are still small. 

When the basis is not complete, the calculated single 
particle energies contain the contributions from the other 
levels in both the Fermi sea and the Dirac sea. The 
contributions from the Fermi sea raises the total energy, 
while that from the Dirac sea lowers the energy. Thus 
as more basis states are included in the expansion, the 
binding energy may not increase monotonically. This 
is different from the non-rclativistic calculations where 
the calculated binding energy always increases when the 
basis size is increased. In the whole deformation range 
considered here, the largest difference between the results 
with A'^ph = 16 and A^ph = 20 appears at the top of 
the second barrier. This difference is about 300 keV. 
It is very time consuming if non-axial deformations are 
included in the calculations. In the present work, for 
the PES's of actinidc nuclei we use A^ph = 20 in axially 
symmetric calculations and A^ph = 16 when the triaxial 
deformations are included. 



basis 



FIG. 2. (Color online) Calculated total binding energy of 
240py by using basis with different basis deformation /3basis 
and different number of major shells A'ph. The results cal- 
culated with 16, 18, and 20 shells are denoted by red, green, 
and blue dots, respectively. The deformation is constrained 
to /?20 = 0.3, 1.3, 2.0, and 3.0, respectively. 



cal two-humped structure and the reflection asymmetric 
shapes are favored in the outer barrier region. To see 
more clearly the truncation errors we have amplified the 
figure near four important points, i.e., the ground state, 
the top of the inner barrier, the isomeric state, and the 
top of the outer barrier. From the four sub-figures we 
can investigate the convergence properties of our model. 
When A^ph increases from 16 to 20, the binding energy 
changes less than 200 keV for the ground state and less 
than 150 keV for the second minimum. The changes of 
the first and the second barrier heights are within 200 
keV and 300 keV, respectively. 

Next we show how the results depend on the deforma- 
tion of the ADHO potential which is used to generate 
the ADHO basis; for brevity, we will call it the basis de- 
formation and label it with /3basis- In Fig- 2 we depict 
the calculated binding energy as a function of the basis 
deformation. In principle, if the basis space is complete 
or nearly complete, the results should not change a lot 
when the basis deformation changes. Near the ground 
state /320 = 0.3, the calculated energies are rather stable 



B. Three-dimensional PES of ^*''Pu 

The double-humped fission barriers of actinide nuclei 
arc usually used as a benchmark for theoretical mod- 
els [91, 96, 98, 99, 112, 113, 153-160]. As a nucleus 
evolves from the ground state to the fission configura- 
tions, various shape degrees of freedom play important 
and different roles in determining the heights of the inner 
and outer barriers. It has long been known that the in- 
ner barrier is lowered when the triaxial deformation is al- 
lowed, while for the outer barrier, the reflection asymmet- 
ric shape is favored. We have shown that the reflection 
asymmetric outer barrier may be further lowered by in- 
cluding the non-axial shape degrees of freedom [125-127]. 
In this section we show some results of multi-dimensional 
PES's for actinides. We will flrst present detailed results 
about the three dimensional PES of -^"^^Pu. then display 
systematic results of the actinide nuclei. In these calcu- 
lations, the effective interaction PC-PKl [152] is used. 

In Figs. 3 and 4 is shown the calculated three dimen- 
sional PES of ^'*"Pu. In each sub-figure we fix the value 
of /320 and display the energy as a function of /322 and /^sq. 
In other words, we have made three-dimensional con- 
straints on the corresponding multipole moments. Note 
that other shape degrees of freedom /3a^ with even /i, 
e.g., /332, /340, P42, PiA, Pso, ■■■, are also included in 
the calculations self-consistently. The MDC-RMF equa- 
tions are solved for each point on the deformation lattice 
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FIG. 3. (Color online) Sections of the three-dimensional potential energy surfaces E = i?(/320,/322,/?3o) of ^''"Pu near the inner 
barrier from MDC-RMF calculations. In each sub-figure the energy is shown as a function of the deformation parameters /322 
and /?3o when /320 is fixed to certain values. The energies are normalized with respect to the binding energy of the ground state. 
The contour interval is 1 MeV. 
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FIG. 4. (Color online) Sections of the three-dimensional potential energy surfaces E — E{l32o, 1322, Pso) of '^'"'Pu near the outer 
barrier from MDC-RMF calculations. In each sub-figure the energy is shown as a function of the deformation parameters /322 
and /J30 when /320 is fixed to certain values. The energies are normalized with respect to the binding energy of the ground state. 
The contour interval is 1 MeV. 



(/?20:/322,/33o) in which /32o runs from 0.25 to 0.95 with 
a step size of 0.05, /322 from to 0.25 with a step size 
of 0.01, and fS-so from to 0.50 with a step size of 0.05. 
The points with f322{P3o) < are obtained through the 
relation £;(^20,^22,/33o) = ^(/320, |/322|, l/SsoD- That is, 
for plotting each sub-figure 26 (for /322) x 11 (for P^q) 
= 286 points arc calculated, and totally 286 x 15 (for 
/32o) = 4290 points are calculated for Figs. 3 and 4. This 
deformation lattice covers the shape space of the most 
interest for ^*°Pu, from the ground state to the isomeric 
state and fissioning configurations. 

Now let us examine the first three sub-figures with 
1^20 = 0.25, 0.30, and 0.35. It is clear that the ground 
state with /?2o '^ 0.3 is both axial symmetric and re- 
flection symmetric, though it is a little soft against the 
octupole distortion. When the nucleus is stretched by 
the quadrupole constraining potential, it becomes softer 
against the triaxial as well as the octupole distortions. 
From the sub-figures with /32o = 0.40 ^ 0.65 one finds 
that two symmetric minima with non-zero triaxial defor- 
mation /322 appear and the corresponding fission paths 
are much more favored than the axial symmetric one, 
which is consistent with earlier calculations [100]. Al- 
though Pso = for all minima in these sub-figures, the 
softness against the octupole distortion changes as the 
nucleus is elongated. The inner fission barrier locates 
near the deformation /32o ^ 0.60. For the last five sub- 
figures with /320 — 0.70 ^ 0.95; the situation becomes 
more simpler, where the nucleus becomes axial symmet- 
ric and reflection symmetric again. Nevertheless, the 
trend to become softer against the octupole distortion 
can be seen from these sub-figures, which indicates that 
the refiection asymmetric shape becomes more relevant. 

As the deformation /320 becomes larger, the second 
minimum and the second saddle point of the PES ap- 
pear. From the last three sub-figures of Fig. 3 and the 
first three sub-figures of Fig. 4 we see that, the nucleus 
keeps refiection symmetric near the fission isomeric state 
with /320 ^ 0.95; but it becomes softer against the /S^q 
distortion. Note that the scales of the Pso axes are dif- 
ferent in Figs. 3 and 4. At /32o = 1.15, two minima corre- 
sponding to reflection asymmetric shapes appear. Here 
the effects of the non-axial deformation is not apparent, 
but along the /322 direction the PES becomes softer. Fi- 
nally, at /320 = 1.2 the energy of the reflection asymmetric 
shape is lower by about 1 MeV than that of the reflec- 
tion symmetric one. Interestingly, when /32o increases 
continuously, beyond the top of the second barrier with 
/320 = 1-2, each reflection asymmetric minimum splits 
into two minima with non-zero /?22- The energy gain due 
to the triaxial distortion is about 1 MeV. The nucleus 
becomes axial symmetric again when /32o > 1-4. 

Prom above discussions, we can draw the following con- 
clusions for ^^°Pu: (1) Both the ground state and the 
fission isomeric state are axial and reflection symmetric; 
(2) Around the first fission barrier it assumes triaxial and 
reflection symmetric shapes; (3) Around the second fis- 
sion barrier both triaxial and octupole deformations are 



important. 



C. PES's of even-even actinide nuclei around 
fission barriers 

The self-consistent three-dimensional constraint calcu- 
lations are very time-consuming, we have only performed 
such calculations for ^""^Pu. From this benchmark study 
we learned many experiences about the important roles 
played by various shape degrees of freedom in different re- 
gions of the deformation space, including the conclusions 
listed in the end of Sec. Ill B. These experiences are used 
in a systematic study of even-even actinide nuclei. 

Since around the inner barrier an actinide nucleus as- 
sumes triaxial but reflection symmetric shapes, we can 
make a one-dimensional constraint calculation with the 
triaxial deformation allowed and the reflection symmetry 
imposed. In Fig. 5 we present potential energy curves of 
even-even actinide nuclei around the flrst fission barrier 
from MDC-RMF calculations. The axial symmetric re- 
sults are displayed by dashed lines and those from triaxial 
calculations are shown by solid lines. The empirical val- 
ues of fission barrier heights are taken from Ref. [161] and 
labelled with solid red circles. We also list in Table I the 
heights of the first and second fission barriers of actinide 
nuclei from MDC-RMF calculations compared with the 
empirical values. It is clearly seen that the triaxial defor- 
mation lowers the inner barrier of these actinide nuclei by 
about 1 '--^ 4 MeV and this is consistent with Ref. [100]. 
The agreement of our calculation results with the empir- 
ical ones is good for most of the nuclei studied here with 
exceptions in the two thorium isotopes and 238,240y^ ^g 
have discussed possible reasons for these discrepancies in 
Ref. [125]. 

It is more complicated to calculate the second fission 
barrier height because more shape degrees of freedom be- 
come important and there are often two or more fission 
paths. What we have done is the following [125]: (1) 
The axial symmetry is assumed and a two-dimensional 
constraint calculation in the {l320iP3o) plane is made; (2) 
From the two-dimensional calculation the lowest fission 
path f^so'"^^ {P20) is approximately identified; (3) Along 
this fission path, a one-dimensional /32o-constraint calcu- 
lation is performed with the triaxial and octupole de- 
formations allowed; (4) At each point with /320 fixed, 
the initial deformations are taken as f3^^^' = and 
/3™- = /3^°o™"'(/32o); (5) In this one-dimensional poten- 
tial energy curve, the second saddle point is located and 
the second barrier height is extracted. 

Potential energy curves of even-even actinide nuclei 
around the second fission barrier from MDC-RMF cal- 
culations are shown in Fig. 6 from which and Table I 
one finds that for most of the nuclei investigated here, 
the triaxiality lowers the second barrier by about 0.5 
~ 1 MeV, accounting for about 10 ~ 20% of the bar- 
rier height. The calculation results with the triaxiality 
included agrees well with the empirical values for all ac- 
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FIG. 5. (Color online) Potential energy curves of even-even actinide nuclei around the first fission barrier from MDC-RMF 
calculations. The axial symmetric results are displayed by dashed lines, while those from the triaxial calculations are shown by 
solid lines. The binding energies are normalized with respect to the binding energies of the ground states of each nucleus. The 
empirical values of fission barriers are taken from Ref. [161]. 
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FIG. 6. (Color online) Potential energy curves of even-even actinide nuclei around the second fission barrier from MDC-RMF 
calculations. The reflection asymmetric shapes are allowed in the calculations. The axial symmetric results are displayed by 
dashed lines, while those from the triaxial calculations are shown by solid lines. The binding energies are normalized with 
respect to the binding energies of the ground states of each nucleus. The empirical values of fission barriers are taken from 
Ref. [161]. 
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TABLE I. Heights of the first and second fission barriers of 
some even-even actinide nuclei from MDC-RMF calculations 
compared with the empirical values ( "Emp" ) which are taken 
from Ref. [161]. The calculation results with and without the 
axial symmetry imposed are denoted by "AS" (axial symmet- 
ric) and "TA" (triaxial), respectively. The height is in MeV. 



TABLE IL Heights of the second barrier of ^^"^U calculated 
from MDC-RMF models with different parameter sets. The 
results with and without non-axial deformations included and 
their difference are denoted by Bas, Bta, and AB, respec- 
tively. 



Nucleus 


Z 


N 


A 


First bai 


■rier 


Second barrier 










AS TA 


Emp 


AS 


TA 


Emp 


^^"Th 


90 


140 


230 


5.03 3.96 


6.10 


6.80 


6.37 


6.50 


232tj^ 


90 


142 


232 


4.94 4.12 


5.80 


6.70 


6.18 


6.70 


232 u 


92 


140 


232 


5.71 4.81 


4.90 


6.20 


5.64 


5.40 


234u 


92 


142 


234 


6.15 5.09 


4.80 


6.20 


5.55 


5.50 


236 u 


92 


144 


236 


6.40 5.11 


5.00 


6.15 


5.31 


5.70 


238 u 


92 


146 


238 


6.54 5.03 


6.30 


6.20 


5.42 


5.50 


240 u 


92 


148 


240 


6.58 4.96 


6.10 


6.38 


5.43 


5.80 


238p^ 


94 


144 


238 


7.72 5.96 


5.60 


6.05 


5.56 


5.10 


240 p^ 


94 


146 


240 


7.98 5.92 


6.05 


6.24 


5.60 


5.20 


242p^ 


94 


148 


242 


8.05 5.77 


5.85 


6.43 


5.74 


5.10 


244 p^ 


94 


150 


244 


7.85 5.40 


5.70 


6.26 


5.49 


4.90 


246p^ 


94 


152 


246 


7.37 4.76 


5.40 


5.84 4.96 


5.30 


2*2Cm 


96 


146 


242 


8.80 6.49 


6.65 


5.72 


4.85 


5.00 


2**Cm 


96 


148 


244 


9.04 6.34 


6.18 


5.90 


4.88 


5.10 


246 Cm 


96 


150 


246 


8.89 5.84 


6.00 


5.40 


4.62 


4.80 


2«Cm 


96 


152 


248 


8.43 5.35 


5.80 


4.10 


— 


4.80 


2S0Cm 


96 


154 


250 


7.77 4.79 


5.40 


2.60 


— 


4.40 


250 (.f 


98 


152 


250 


8.87 5.70 


5.60 


2.40 


— 


3.80 


252 cf 


98 


154 


252 


8.41 5.26 


5.30 


1.20 


— 


3.50 



tinidc nuclei shown in Fig. 6. In Ref. [125], we have found 
that for ^"'^Cni, the height of the second barrier without 
the triaxial deformation included is already smaller than 
the empirical value. From Table I it is seen that this 
is also the case for ^^°Cm and 250,252q£ Therefore the 
heights of the second barriers with the triaxial deforma- 
tion included are not listed for these nuclei. The reason 
for these discrepancies may be related to a strong com- 
petition between the two or among more fission paths 
and the assumption on the barrier width made when the 
empirical value of the barrier height is evaluated [125]. 

From Figs. 5 and 6 one can find some interesting fea- 
tures concerning the triaxial shape of actinides around 
the first and the second barriers. Though it is emphasized 
here and in Refs. [100, 101] that the triaxiality is crucial 
in determining the height of these barriers, it it not nec- 
essary that the nucleus in question is triaxially deformed 
at these two saddle points. One can find that all of these 
even-even actinide nuclei shown in Fig. 5 are axially de- 
formed at the first saddle points from our MDC-RMF cal- 
culations. From Fig. 6 it is seen that except 232,234,240jj^ 
all other nuclei are also axially deformed at the second 
saddle points. The triaxial shape is only important be- 
fore the saddle point around the first barriers, while it is 
often so after the saddle point around the second barrier, 
with exceptions in uranium isotopes. Therefore the two 
saddle points become closer when the triaxial deforma- 
tion is included. Further investigations should be made 
to reveal the physics behind these features. 



Parameter set 


Bas (MeV) Bta (MeV) 


AB (MeV) 


NL3* 


7.54 


6.85 


0.69 


NL-Z2 


4.83 


3.91 


0.92 


PC-PKl 


6.20 


5.55 


0.65 


DD-ME2 


8.19 


7.51 


0.68 


DD-PCl 


6.13 


5.64 


0.49 



D. Parameter dependence of the effects of 
triaxiality on the second barrier height 

From the above detailed studies of the PES's of ac- 
tinide nuclei in the (/320i f^22, /^so) space, it is clear that 
the triaxiality plays an important role upon the second 
fission barriers of actinide nuclei. We have studied the 
parameter dependence of the effects of triaxiality on the 
second barrier height and some results can be found in 
Refs. [125, 126, 162]. In Table II are listed the second 
barrier heights of ^'^^U calculated with different parame- 
ter sets, including meson-exchange ones NL3* [163], NL- 
Z2 [164], DD-ME2 [165], and point-coupling ones PC- 
PKl [152] and DD-PCl [166]. One finds that the sec- 
ond barrier height may differ by a few MeV with dif- 
ferent effective interactions, but in all cases the barriers 
are lowered considerably by the non-axial deformations. 
For this specific nucleus, the lowering effects with differ- 
ent effective interactions are about at least 0.5 MeV and 
may as large as 0.9 MeV which are clearly larger than 
the possible errors from the basis truncations discussed 
in Sec. Ill A. 



E. Further discussions 

A systematic study of fission barriers is performed for 
even-even superheavy nuclei with Z = 112 ~ 120 by 
using covariant density functional models and the outer 
fission barriers are found to be considerably affected by 
the triaxiality [101]. However, in a recent work within the 
macroscopic-microscopic model, it was found that the in- 
fluence of the triaxiality on second fission barriers is neg- 
ligible [167]. This raises an open problem: What are the 
main reasons for the different effects of the triaxiality on 
the second fission barriers predicted by the MM model 
and CDFTs? It would be interesting to make more in- 
vestigations with different models, especially, with the 
non-relativistic DFTs. 

One of the problems concerning the PES calculated 
from self-consistent approaches is that there may exist 
unexpected discontinuities on the PES's [73]. This is 
mainly due to the complexity of multi-dimensional PES's. 
When a high-dimensional PES is projected onto a low- 
dimensional one, the minimizing procedure incorporated 
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implicitly in the self-consistent approaches may cause a 
sudden change in the total energy at some points in the 
low-dimensional deformation space. This may result in 
spurious saddle points. In our MDC-RMF calculations, 
we have tried some ways to avoid discontinuities on the 
PES's and to exclude spurious saddle points [125]. To 
be completely free of discontinuities or spurious saddle 
points, one certainly should carry out multi-dimensional 
constraint calculations with higher-polarity deformations 
included. In Ref. [168] the authors analyzed the origin of 
the discontinuities and proposed a numerical method to 
identify them. It will be useful to implement this method 
in the present MDC-RMF calculations; this will be one 
of our topics in the future. 

Since many symmetries are broken in mean field cal- 
culations, quantum numbers related to these symmetries 
are not conserved any more [2]. For example, the to- 
tal angular momentum or nuclear spin J is not a good 
quantum number when /32o is not zero; similarly, the pro- 
jection of the total angular momentum on the symme- 
try axis K is also not conserved if a nucleus is triaxi- 
ally deformed [169]. When the reflection symmetry is 
broken, the parity is not a good quantum number. In 
our three-dimensional constraint calculations, all these 
quantum numbers arc not good ones. One may intro- 
duce techniques of projection to restore the broken sym- 
metries [2]. In recent years, the importance of various 
projections on the PES's and fission barriers have been 
investigated [94, 160, 170, 171]. 



IV. 



SUMMARY 



We developed multi-dimensional constraint rclativistic 
mean field (MDC-RMF) models. In these models, the nu- 
clear shape is assumed to be invariant under the reversion 
of X and y axes, i.e., the intrinsic symmetry group is V4 
and all shape degrees of freedom j3x^ with even /i (/320i 
/322, /330, /332, /340, ' ' ' ) are included self-consistently. The 
RMF functional can be one of the following four forms: 
the meson exchange or point-coupling nuclcon interac- 
tions combined with the non-linear or density-dependent 
couplings. We solve the Dirac equation in an axially de- 
formed harmonic oscillator (ADHO) basis. The conver- 
gence of the calculated results against the basis trunca- 
tion and the basis deformation is studied and it is shown 
that reasonably large ADHO basis is able to provide a 
desired accuracy. 

Three-dimensional potential energy surface in the (/320, 
/322, /330) plane is obtained for ^"'"Pu and potential energy 
curves of actinide nuclei around the first and second fis- 
sion barriers are studied systematically. It is found that 
the triaxiality is crucial in determining the height of both 
the first and the second fission barriers. However, most 
of actinides studied here are axially deformed at the first 
and second saddle points. The triaxial shape is mostly 
important before the saddle point around the first bar- 
rier, while it is often so after the saddle point around the 



second barrier. Taking ^'^^U as an example, we have stud- 
ied the parameter dependence of the effects of triaxiality 
on the second barrier height and found that the height 
of the second barrier may differ by a few MeV with dif- 
ferent effective interactions, but in all cases the barriers 
are lowered considerably by the non-axial deformations. 
We conclude that it is important to include the reflec- 
tion asymmetric and non-axial shapes simultaneously for 
the study of potential energy surfaces and fission barriers 
of actinide nuclei and of those in unknown mass regions 
such as superheavy nuclei. 
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Appendix A: Matrix elements of the Dirac 
Hamiltonian 



In the ADHO basis, the Dirac equation (5) is trans- 
formed into a matrix eigenvalue problem, 

5](a|M+(r)|a')/r' +EH'^-pI«')5?' -eJf, 

a' a' 

^(a|<T . p\a')ff + Y.{a\M^{r)\a')gf = e,gf . 

a' a' 

In the cylindrical coordinate system the kinetic energy 
operator a ■ p reads. 



CT ■ p = I 






e-'v{dp-ip-^d^) 



-d. 



The matrix elements of the kinetic energy operator a ■ p 
can be calculated as: 

{a\(T -pla') ^ {n^,np,mi,ms\(T ■ p\n'^,n'p,m'i,m'^) 

+ (5„^ _i(5„'^,i(5„,,«i {Baa' - m'lCaa') 
+ ^m„^^m'^,-^^n,,n'^ {Baa' + m'lCaa') , 

(Al) 
where the integrals read, 

Aaa' = / dz4>n,dz4>n', (A2) 
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(A3) 
(A4) 



p—f,gaa' \i—l / 

Note that the factors with ip in $„ and $„/ arc in an 
exponential form, rearranging the terms, pvif) can be 
written as, 

U{r)^U{z,p,if) = — Y^ U^^'\z,p)e'^'^, (A5) Pv(r) = — U^ (^,p) + 2 ^ ^^^"^(z, p) cos(2nv?) , 

(B2) 



The matrix elements of the potentials U{r) — M±{r) 
are calculated by expanding U{r) into a Fourier series. 



where the components U^^^Zjp) are calculated by the 
inverse Fourier transformation, 

U^^\z,p) = f "d^C/(z,p,^)e-'^^, (A6) 



where the components 



P'^' 



E E [f2'^h?pA5K'-K,,s^ 



where the integral over ip is performed on a uniformly 
distributed mesh points. From the symmetry conditions 
it is easy to find, 

[/(m) ^[/(-m) = (_i)Mf7(-M). (A7) 

Thus the Fourier series is abbreviated as 



xc:c„.0„>„.i?:r;C 



(B3) 



arc calculated with the wave functions /" and gf. 

The scalar density psir) can be calculated similarly. 
The isovector densities are just the linear combinations 
of the corresponding densities of protons and neutrons. 

The derivatives of the vector density reads 



Uir) = ^(u^'Hz,p) + 2f^U^^-\z,p)cos{2ncp)y V'pv{r)^ E E ( E^''^"?'" j ^'(*"*"')' (^4) 

\ n=l / P=f,g aa' \i=l / 

(A8) 



Finally the matrix element of U{r) reads. 



{a\U\a') = {nz,nr,mi,ms\U\n'^,n'^,m'i,m'^] 

1 r, 

2^ 









^ {5K'-K,2n + 5k-K', 2n) -D^^, 



p=f,g ' 

where the derivatives of the basis can be calculated as, 

V2($t^$„,) = /„„, + 2J„„, + /:,„, (B5) 

with laa' and Jaa' defined later. First, from the eigen 
equation Eq. (10) we get 



CA9) 



where we have performed the integrals over if and 



thus 



oo />oo 



dz / pdpf/('^)(z,p)</)„,</)„-i?™'i?:':'. 



J-oo JO 

(AlO) 
are performed numerically by using the Gaussian quadra- 
ture. For axial symmetric potentials only the first term 
survives, while for triaxially deformed nucleus we must 
calculate additional terms with n ^ 0. 



Appendix B: Densities and their derivatives 

After solving the Dirac equation, we can calculate the 
densities from /" and g". All the densities are lin- 
ear combinations of the quantities pf{r,T) and pg{r,T) 
which are the contributions to the vector density from 
the large and small components, respectively, r repre- 
sents the isospin. For vector potential we have 

N 

«=i p=f,a 



V2$„ = -2A/ K - Vb{z,p)] $a, (B6) 



= -2M [E^ - Vb{z, p)] (5„,,„,. C:a. 

X (/)„>„' RTR7 X -^e^^™:-™')'^. (B7) 
Second, 



0„.0„: [dpRZypR';:l+mim[p- 



2 f>mi D™i 



KlRu'dAn^dAu'^ 



X —e'C™;-™')"^. (B8) 
27r 



Substituting them into Eq. (B5) and combining the same 
exponentials of ip, we get 

VV(r) = i- fpf (z,p) + 2f;p(,^")(z,p)cos(2n^)j , 

(B9) 
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in which 



p—f,g OLO.' \i— 1 / 



(BIO) 
are the Fourier components and 



D^^, (z, p) = ^^^ - 2M [Ea,, +E^~ 2Vb{.z, p)] . 
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